#include <iostream>
#include <cmath>

class EquationSolver
{
public:
  virtual double solve(double a, double b) = 0;
};

class BisectionSolver : public EquationSolver
{
public:
  BisectionSolver(double (*f_)(double)) : f(f_)
  {
    epsilon = 1e-6;
    delta = 1e-6;
    M = 10000;
  };
  BisectionSolver(double (*f_)(double), double Epsilon, double Delta, int m) : f(f_), epsilon(Epsilon), delta(Delta), M(m) {};
  
  double solve(double a, double b) override
  {
    int k;
    double c, w;
    double h = b - a, u = f(a);
    for( k = 0; k <= M; k++)
      {
	h = h / 2;
	c = a + h;
	if( h < delta || k == M)
	  break;
	w = f(c);
	if( w < epsilon && w > -epsilon )
	  break;
	else if( w*u > 0 )
	  a = c;
      }
    return c;
  };
private:
  double (*f)(double);
  double epsilon, delta;
  int M;
};

class NewtonSolver : public EquationSolver
{
public:
  NewtonSolver(double (*f_)(double), double (*fp_)(double)) : f(f_), fp(fp_)
  {
    epsilon = 1e-6;
    delta = 1e-6;
    M = 10000;
  };
  NewtonSolver(double (*f_)(double), double (*fp_)(double), double Epsilon, double Delta, int m) : f(f_), fp(fp_), epsilon(Epsilon), delta(Delta), M(m) {};
  
  double solve(double a, double b) override
  {
    double x = (a + b) / 2, u;
    int k;
    for( k = 0; k <= M; k++ )
      {
	u = f(x);
	if( u < epsilon && u > -epsilon )
	  break;
	x = x - u / fp(x);
      }
    return x;
  };
private:
  double (*f)(double);
  double (*fp)(double);
  double epsilon, delta;
  int M;
};

class SecantSolver : public EquationSolver
{
public:
  SecantSolver(double (*f_)(double)) : f(f_)
  {
    epsilon = 1e-10;
    delta = 1e-10;
    M = 10000;
  };
  SecantSolver(double (*f_)(double), double Epsilon, double Delta, int m) : f(f_), epsilon(Epsilon), delta(Delta), M(m) {};
  
  double solve(double x0, double x1) override
  {
    double a = x0, b = x1, s, h;
    double u = f(b), v = f(a);
    int k;
    for( k = 2; k <=M; k++ )
      {
	s = (b - a) / (u - v);
	a = b;
	v = u;
	b = b - u * s;
	h = b - a;
	if( h < delta && h > -delta )
	  break;
	u = f(b);
	if( u < epsilon && u > -epsilon )
	  break;
      }
    return b;
  };
private:
double (*f)(double);
  double epsilon, delta;
  int M;
};

double fB_1(double x)
{
  return 1 / x - tan(x);
};

double fB_2(double x)
{
  return 1 / x - pow(2,x);
};

double fB_3(double x)
{
  return pow(2,-x) + exp(x) + 2 * cos(x) - 6;
};

double fB_4(double x)
{
  return (pow(x,3) + 4 * pow(x,2) + 3 * x + 5) / (2 * pow(x,3) - 9 * pow(x,2) + 18 * x - 2);
};

double fC(double x)
{
  return x - tan(x);
};

double fpC(double x)
{
  return -1 / pow(cos(x),2);
};

double fD_1(double x)
{
  return sin(x/2) - 1;
};

double fD_2(double x)
{
  return exp(x) - tan(x);
};

double fD_3(double x)
{
  return pow(x,3) - 12 * pow(x,2) + 3 * x + 1;
};

double fE(double x)
{
  double V = 12.4, L = 10, r = 1;
  return V / L + r * r * asin(x/r) + x * sqrt(r*r-x*x) - 0.5 * M_PI * r * r; 
};

double fpE(double x)
{
  double V = 12.4, L = 10, r = 1;
  return 2 * sqrt(r*r-x*x); 
};

double fF_a(double alpha)
{
  double beta = M_PI*11.5/180, l = 89, h = 49, D = 55;
  double A = l*sin(beta), B = l*cos(beta), C = (h+D/2)*sin(beta)-D*tan(beta)/2, E = (h+D/2)*cos(beta)-D/2;
  return A*sin(alpha)*cos(alpha)+B*sin(alpha)*sin(alpha)-C*cos(alpha)-E*sin(alpha); 
}

double fpF_a(double alpha)
{
  double beta = M_PI*11.5/180, l = 89, h = 49, D = 55;
  double A = l*sin(beta), B = l*cos(beta), C = (h+D/2)*sin(beta)-D*tan(beta)/2, E = (h+D/2)*cos(beta)-D/2;
  return A*cos(2*alpha)+2*B*sin(alpha)*cos(alpha)+C*sin(alpha)-E*cos(alpha); 
}

double fF_b(double alpha)
{
  double beta = M_PI*11.5/180, l = 89, h = 49, D = 30;
  double A = l*sin(beta), B = l*cos(beta), C = (h+D/2)*sin(beta)-D*tan(beta)/2, E = (h+D/2)*cos(beta)-D/2;
  return A*sin(alpha)*cos(alpha)+B*sin(alpha)*sin(alpha)-C*cos(alpha)-E*sin(alpha); 
}

double fpF_b(double alpha)
{
  double beta = M_PI*11.5/180, l = 89, h = 49, D = 30;
  double A = l*sin(beta), B = l*cos(beta), C = (h+D/2)*sin(beta)-D*tan(beta)/2, E = (h+D/2)*cos(beta)-D/2;
  return A*cos(2*alpha)+2*B*sin(alpha)*cos(alpha)+C*sin(alpha)-E*cos(alpha); 
}
  
int main()
{
  EquationSolver* solverB_1;
  solverB_1 = new BisectionSolver(fB_1);
  double bisectionResultB_1;
  bisectionResultB_1 = solverB_1->solve(0.0,M_PI/2);
  std::cout << "Problem B: " << std::endl;
  std::cout << "1. " << "alpha = " << bisectionResultB_1 << ", f(alpha) = " << fB_1(bisectionResultB_1) << std::endl;

  EquationSolver* solverB_2;
  solverB_2 = new BisectionSolver(fB_2);
  double bisectionResultB_2;
  bisectionResultB_2 = solverB_2->solve(0.0,1.0);
  std::cout << "2. " << "alpha = " << bisectionResultB_2 << ", f(alpha) = " << fB_2(bisectionResultB_2) << std::endl;

  EquationSolver* solverB_3;
  solverB_3 = new BisectionSolver(fB_3);
  double bisectionResultB_3;
  bisectionResultB_3 = solverB_3->solve(1.0,3.0);
  std::cout << "3. " << "alpha = " << bisectionResultB_3 << ", f(alpha) = " << fB_3(bisectionResultB_3) << std::endl;

  EquationSolver* solverB_4;
  solverB_4 = new BisectionSolver(fB_4);
  double bisectionResultB_4;
  bisectionResultB_4 = solverB_4->solve(0.0,4.0);
  std::cout << "4. " << "alpha = " << bisectionResultB_4 << ", f(alpha) = " << fB_4(bisectionResultB_4) << std::endl;

  EquationSolver* solverC;
  solverC = new NewtonSolver(fC,fpC);
  double newtonResultC_1, newtonResultC_2;
  std::cout << "Problem C: " << std::endl;
  newtonResultC_1 = solverC->solve(4.0,5.0);
  std::cout << "the alpha near 4.5, alpha =  " << newtonResultC_1 << "f(alpha) = " << fC(newtonResultC_1) << std::endl;
  newtonResultC_2 = solverC->solve(7.6,7.8);
  std::cout << "the alpha near 7.7, alpha =  " << newtonResultC_2 << "f(alpha) = " << fC(newtonResultC_2) << std::endl;

  SecantSolver* solverD_1;
  solverD_1 = new SecantSolver(fD_1);
  double secantResultD_1;
  secantResultD_1 = solverD_1->solve(0.0,M_PI/2);
  std::cout << "Problem D: " << std::endl;
  std::cout << "1. " << "alpha = " << secantResultD_1 << ", f(alpha) = " << fD_1(secantResultD_1) << std::endl;

  SecantSolver* solverD_2;
  solverD_2 = new SecantSolver(fD_2);
  double secantResultD_2;
  secantResultD_2 = solverD_2->solve(1.0,1.4);
  std::cout << "2. " << "alpha = " << secantResultD_2 << ", f(alpha) = " << fD_2(secantResultD_2) << std::endl;

  SecantSolver* solverD_3;
  solverD_3 = new SecantSolver(fD_3);
  double secantResultD_3;
  secantResultD_3 = solverD_3->solve(0.0,-0.5);
  std::cout << "3. " << "alpha = " << secantResultD_3 << ", f(alpha) = " << fD_3(secantResultD_3) << std::endl;

  EquationSolver* solverE_1;
  solverE_1 = new BisectionSolver(fE);
  double bisectionResultE_1;
  bisectionResultE_1 = solverE_1->solve(0.0,1.0);
  std::cout << "Problem E: " << std::endl;
  std::cout << "1. " << "h = " << bisectionResultE_1 << ", f(h) = " << fE(bisectionResultE_1) << std::endl;

  EquationSolver* solverE_2;
  solverE_2 = new NewtonSolver(fE,fpE);
  double newtonResultE_2;
  newtonResultE_2 = solverE_2->solve(0.0,1.0);
  std::cout << "2. " << "h = " << newtonResultE_2 << ", f(h) = " << fE(newtonResultE_2) << std::endl;

  SecantSolver* solverE_3;
  solverE_3 = new SecantSolver(fE);
  double secantResultE_3;
  secantResultE_3 = solverE_3->solve(0.0,1.0);
  std::cout << "3. " << "h = " << secantResultE_3 << ", f(h) = " << fE(secantResultE_3) << std::endl;

  EquationSolver* solverF_a;
  solverF_a = new NewtonSolver(fF_a,fpF_a);
  double newtonResultF_a;
  std::cout << "Problem F: " << std::endl;
  newtonResultF_a = solverF_a->solve(0.0,1.0)*180/M_PI;
  std::cout << "a. " << "alpha = " << newtonResultF_a << ", f(alpha) = " << fF_a(newtonResultF_a*M_PI/180) << std::endl;

  EquationSolver* solverF_b;
  solverF_b = new NewtonSolver(fF_b,fpF_b);
  double newtonResultF_b;
  newtonResultF_b = solverF_b->solve(0.0,1.0)*180/M_PI;
  std::cout << "b. " << "alpha = " << newtonResultF_b << ", f(alpha) = " << fF_b(newtonResultF_b*M_PI/180) << std::endl;

  SecantSolver* solverF_c;
  solverF_c = new SecantSolver(fF_a);
  double secantResultF_c;
  secantResultF_c = solverF_c->solve(1.0,3.0)*180/M_PI;
  std::cout << "c. " << "alpha = " << secantResultF_c << ", f(alpha) = " << fF_a(secantResultF_c*M_PI/180) << std::endl;

  return 0;
}
